A transcriptome-informed QSP model of metastatic triple-negative breast cancer identifies predictive biomarkers for PD-1 inhibition

Triple-negative breast cancer (TNBC), a highly metastatic breast cancer subtype, has limited treatment options. While a small number of patients attain clinical benefit with single-agent checkpoint inhibitors, identifying these patients before the therapy remains challenging. Here, we developed a transcriptome-informed quantitative systems pharmacology model of metastatic TNBC by integrating heterogenous metastatic tumors. In silico clinical trial with an anti–PD-1 drug, pembrolizumab, predicted that several features, such as the density of antigen-presenting cells, the fraction of cytotoxic T cells in lymph nodes, and the richness of cancer clones in tumors, could serve individually as biomarkers but had a higher predictive power as combinations of two biomarkers. We showed that PD-1 inhibition neither consistently enhanced all antitumorigenic factors nor suppressed all protumorigenic factors but ultimately reduced the tumor carrying capacity. Collectively, our predictions suggest several candidate biomarkers that might effectively predict the response to pembrolizumab monotherapy and potential therapeutic targets to develop treatment strategies for metastatic TNBC.


INTRODUCTION
Numerous hallmarks of cancer that facilitate the transformation of normal cells to neoplastic cells have been proposed (1). Metastasis, a process whereby cancer cells acquire invasive ability and spread from the site of origin to distant organs, is one of the hallmarks and accounts for 90% of cancer-related deaths (2,3). While the metastatic potential varies for different cancer types (4,5), a subtype of breast cancer called triple-negative breast cancer (TNBC) is characterized by a high tendency for metastasis. TNBC constitutes 15 to 20% of breast cancer and is considered a highly aggressive breast cancer subtype (6). TNBC tumors lack expression of estrogen and progesterone receptors, do not overexpress human epidermal growth factor receptor-2 (HER-2), and cannot be treated with drugs that target these pathways in other subtypes of breast cancer. Historically, in terms of systemic therapy, this subtype of breast cancer has relied mostly on chemotherapy in both the early-stage and metastatic settings.
TNBC tumors are usually sensitive to chemotherapy, but many patients will still relapse after treatment and median overall survival is poor (7)(8)(9). For TNBC, the first distant site of metastatic appearance is the lung (40%), followed by brain (30%), liver (20%), and bone (10%) (10). Genomic analysis of metastatic TNBC tumor samples from a rapid autopsy patient suggests that lung metastases promote subclonal diversification that gives rise to highly aggressive cancer clones, thus acting as "subclone incubators" (11).
Reduced immune cell abundance and tumor-infiltrating lymphocytes (TILs) were observed in metastatic tumors compared to primary tumors in breast cancer (12)(13)(14). Discordant programmed cell death-ligand 1 (PD-L1) expression status (12,15) and genomic alterations (16) between metastatic and primary tumors have also been observed. Multiple metastatic tumors in different organs or within the same organ also show substantial heterogeneity in clonal composition (11,17). Lower immunogenicity and intraindividual heterogeneity of metastatic tumors exacerbate the challenges of treating metastatic TNBC.
One of the emerging promising treatment strategies aims to enhance antitumor immune responses by targeting immune checkpoint molecules, such as cytotoxic T lymphocyte-associated antigen 4 (CTLA-4), programmed cell death protein 1 (PD-1), and PD-L1 (18,19). Pembrolizumab, a PD-1 inhibitor, boosts T cell responses by blocking the interactions of PD-1 receptor with its ligands PD-L1 and PD-L2 and has been tested in a number of clinical trials of early-stage and metastatic TNBC (20)(21)(22)(23)(24)(25). When pembrolizumab was combined with chemotherapy in the KEYNOTE-355 trial, progression-free and overall survival were improved in patients with tumors with a PD-L1 combined positive score (CPS) of 10 or greater (25,26), subsequently leading to its approval in the treatment of metastatic PD-L1 CPS ≥ 10 TNBC.
When tested as a single agent in patients with metastatic TNBC, response rates were in the range of 5 to 21% depending on the tumor PD-L1 status and prior treatments received by patients (20)(21)(22)(23). In randomized trials, pembrolizumab monotherapy has shown low response rates, and overall survival was not significantly enhanced compared to chemotherapy (22).
With pembrolizumab monotherapy, objective response rates were lower in previously treated patients compared to untreated patients (20,23), and a subset of patients demonstrated complete and durable response (22). Furthermore, with increasing levels of PD-L1 expression in the tumor, a trend toward increasing response rates and overall survival was observed (22,27). These findings highlight the need to identify predictive biomarkers for clinical benefit from pembrolizumab monotherapy. An extensive analysis of biomarkers in clinical studies presents multiple technical challenges that could be circumvented by exploratory analysis in silico followed by validation in clinical studies.
Quantitative systems pharmacology (QSP) models integrate detailed mechanisms of biological systems with drug pharmacology and assist in drug discovery by promoting a system level understanding of treatment responses (28). QSP models have been developed to investigate the action of various immune checkpoint inhibitors and T cell engagers (29,30), in cancer types such as colorectal cancer (31,32), breast cancer (33)(34)(35)(36), prostate cancer (37), non-small cell lung cancer (38), and melanoma (39). Such models have been used in investigating treatment responses as well as dose optimization and scheduling, have been used in identifying and evaluating biomarkers and preclinical to clinical translation, and were able to explain clinically observed interindividual heterogeneity in treatment responses. As intraindividual heterogeneity is also critical in determining treatment responses in the metastatic setting, Kumar et al. (40) developed a QSP model of metastatic melanoma where multiple heterogenous tumors were explicitly represented, parameterized the model using clinical data from ipilimumab and pembrolizumab monotherapies, and then predicted the response to combination therapy. Model predictions suggested that only tumors with intermediate levels of TILs, referred to as "warm tumors," showed a substantial enhancement in response to combination therapy compared to monotherapies (40).
Because of the high incidence of metastasis in TNBC, understanding the role of intraindividual heterogeneity in treatment responses is critical. In this study, we developed a transcriptomeinformed QSP model with a special focus on lung metastases to investigate patient responses to pembrolizumab monotherapy as a second-or third-line therapy in metastatic TNBC. By performing virtual clinical trials, we characterized the effects of pembrolizumab on the tumor microenvironment in responders and nonresponders and predicted biomarkers potentially associated with clinical benefit.

QSP model development and calibration
We developed a QSP model with metastatic tumors to explore the effects of pembrolizumab monotherapy in metastatic TNBC (Fig. 1). This QSP model uses components of a previously developed QSP model with single TNBC tumor and was extended in this study to incorporate multiple metastatic tumors. Because lung metastases are frequently observed in patients with TNBC and have a unique tumor microenvironment (figs. S1 and S2), we explicitly incorporated lung metastatic tumors in the QSP model by adding two lung tumor compartments that drain into a lung-draining lymph node compartment. In addition, we also incorporated an "other" tumor compartment representing any metastatic tumor other than the lung and an associated other tumor-draining lymph node compartment ( Fig. 1). We calibrated the relative abundance of immune cells in the lung metastatic tumors (with respect to matched primary tumors) in the QSP model using transcriptome-derived relative abundance estimates from lung metastases ( Fig. 2A). Only the abundance estimates of EPIC and quanTIseq were used for model calibration because the estimates are expressed as cell fractions and enable direct comparison with cell fractions from simulations. Similarly, we used the transcriptome-derived relative abundance of cell types in metastatic tumors, excluding the lung, brain, and lymph node samples, to calibrate the immune cell abundance of other metastatic tumors in the QSP model (Fig. 2B). As the availability of transcriptomic data was limited to a low number of tumor samples, the estimates obtained using EPIC and quanTIseq were only treated as rough estimates, and an exact quantitative comparison with simulation results was not performed. We further calibrated the QSP model using data from the KEYNOTE-119 clinical trial (22,27). We generated 1000 virtual patients by selectively varying parameters of the QSP model. Among the 1000 virtual patients, 816 virtual patients who met the target tumor diameter criteria were selected for treatment simulations and subsequent analyses. We assumed that most patients underwent surgery to remove the primary tumor and incorporated this assumption in the model by setting the initial number of cancer cells in the primary tumor compartment to 0. Metastatic tumor compartments were seeded at randomly sampled time points until at least one metastatic tumor reached the target tumor diameter, resulting in a varying number of established metastatic tumors in virtual patients with a maximum of two lung metastatic tumors and one other metastatic tumor. We simulated the action of pembrolizumab monotherapy in virtual patients by administering 200 mg of pembrolizumab every 3 weeks and calibrated the QSP model parameters to recapitulate clinically observed response rates, duration of response, time to response, and the percentage of patients with lung metastases in the virtual clinical trial (Fig. 3A). While a subset of responders showed cancer progression during the course of the treatment, there were also virtual patients who did not progress until the end of the simulation (Fig. 3A, red dots), mimicking patients with long-term response in the clinical trial. As expected, a large intraindividual and interindividual heterogeneity in treatment response was observed among virtual patients (Fig. 3, B to D).

Impact of pembrolizumab treatment on the tumor microenvironment
It is well-known that pembrolizumab enhances the effector function of T cells by blocking the interaction of the PD-1 receptor with its ligands PD-L1 and PD-L2, but the implications of this blockade on the composition and properties of tumor microenvironment are not clear. To characterize this, we compared simulations of virtual patients in the presence and absence of pembrolizumab therapy (fig .  S3). First, we observed an increase in the density of T cells in the central compartment ( fig. S3A), which was expected because of the increased activation of T cells induced by an increase in the cancer cell killing upon the action of pembrolizumab. In the central compartment, an increase in the diversity of cytotoxic T cells and an increase in the ratio of regulatory T cells (T regs ) and cytotoxic T cells were also observed in pembrolizumab-treated virtual patients ( fig. S3, B and C).
Upon pembrolizumab treatment, there was a notable increase in the fraction of immune cells in metastatic tumors including both protumorigenic and antitumorigenic immune cells ( fig. S3E).  Mechanisms considered in the metastatic tumors and associated draining lymph node compartments are the same as mechanisms shown in the primary tumor and associated lymph node compartments, respectively. Th, helper T cell; T reg , regulatory T cell; T cyt , cytotoxic T cell; nTCD4/8, naïve CD4/CD8 T cell; aTCD4/8, activated CD4/8 T cell; IV, intravenous; APC, antigen-presenting cell; mAPC, mature antigen-presenting cell; Mac, macrophage; MDSC, myeloid-derived suppressor cell; Arg-I, arginase I; NO, nitric oxide; MHC, major histocompatibility complex; dLNs, draining lymph nodes.
We then compared the effects of pembrolizumab on immune cell composition and cytokine levels between responders and nonresponders by stratifying virtual patients based on the response status ( Fig. 4 and fig. S5). In general, the aforementioned effects on cellular and molecular species observed in the entire cohort of virtual patients were quantitatively enhanced in responders compared to nonresponders. Increases in the density and diversity of T cells in the central compartment and the fraction of immune cells in the tumor were higher in responders compared to nonresponders ( Fig. 4, A, B, and D). Similarly, the decrease in the richness of cancer clones was greater in responders, suggesting increased clonal extinction (Fig. 4C). A transient decrease in the tumor-infiltrating cytotoxic T cell diversity was also greater in responders (Fig. 4F). Apart from these prominent quantitative differences, mild qualitative differences were observed between responders and nonresponders. While the ratio of protumorigenic and antitumorigenic immune cell densities consistently decreased in nonresponders, a transient increase was observed in responders following an initial decrease (Fig. 4E). Similarly, a transient increase in the density of M1 macrophages and IL-12 levels was seen in responders but was not observed in nonresponders (Fig. 4, G and H). Although IL-2 levels were increased upon pembrolizumab therapy in both responders and nonresponders, the increase tended to be earlier in nonresponders (Fig. 4I).

Identification and ranking of potential biomarkers
We then sought to identify biomarkers that could potentially identify responders to pembrolizumab monotherapy in metastatic TNBC. First, we tested whether the baseline levels of previously identified biomarkers such as tumor PD-L1 expression, the number of neo-antigen-specific T cell clones, and the fraction of TILs correlated with response to pembrolizumab monotherapy. In this analysis, neo-antigen-specific T cell clones were considered as a proxy to tumor mutational burden due to the correlation between these two quantities (41). On the basis of the levels of each of the three biomarkers individually, we generated virtual patient subgroups (see Materials and Methods). We then calculated the fraction of responders among all the patients in each subgroup, hereafter termed as response probability. As expected, we observed a modest increase in response probability among virtual patient subgroups as the threshold levels of these biomarkers were increased (Fig. 5, A to C). Correlation between the PD-L1 expression and fraction of TILs was predicted to be 0.56 ( fig. S6), similar to the clinically reported correlation of 0.45 (42).
We performed an exhaustive testing of biomarker candidates (see Materials and Methods) in the QSP model to search for other biomarkers with a potentially higher predictive power. For each biomarker candidate, we generated virtual patient subgroups and calculated the response probability for each subgroup as previously described. We then ranked the biomarker candidates based on the highest response probability achieved among the virtual patient subgroups generated for each biomarker candidate.
Among the biomarker candidates tested, the density of antigenpresenting cells (APCs) in tumor-draining lymph nodes achieved the highest response probability of 0.86 (Fig. 5D). This was followed by the richness of cancer clones, the fraction of cytotoxic T cells among all T cells in lymph node and central compartments, the fraction of exhausted cytotoxic T cells in metastatic tumors, the tumor diameter, and the diversity and evenness of cytotoxic T cells in lymph nodes. Each of these biomarker candidates achieved a response probability higher than 0.6. Two biomarkers with known clinical relevance, such as tumor PD-L1 expression and the number of neo-antigen-specific T cell clones, were among the best 30 biomarker candidates identified, as both achieved a response probability of approximately 0.5. The fraction of TILs achieved a response probability of 0.45, suggesting that it is a slightly weaker biomarker compared to PD-L1 expression in the tumor.
We tested whether the virtual patient subgroups selected based on identified biomarkers incorporated a considerable number of responders from the whole virtual patient cohort. The virtual patient subgroup chosen based on the density of APCs in lymph nodes, which is the best biomarker identified, included only 11% of responders from the entire virtual patient cohort despite achieving a remarkable response probability (fraction of responders within a selected patient subset rather than the entire patient cohort) of 0.86 ( fig. S7A). In other words, while this biomarker effectively identifies a small subset of responders from the whole cohort (11%), the majority of responders from the whole cohort (89%) are ignored and are not benefitted by this biomarker. Among the virtual patient subgroups chosen based on best 30 biomarkers identified, three biomarkers-the density of T regs in lymph nodes, the fraction of helper and cytotoxic T cells in the central compartment, and the fraction of T regs in the central compartment-also had higher percentage of responders from the whole cohort ( fig. S7A).
To maximize the number of responders in the best subgroups selected, we introduced a metric that we call responder inclusion score (RIS). RIS is calculated as the difference between the fraction of responders and the fraction of nonresponders from the entire cohort present in the subgroup. The first term maximizes the number of responders while the second term minimizes the number of nonresponders in the subgroup with respect to the total number of responders or nonresponders in the entire cohort. When biomarker candidates were ranked on the basis of the RIS rather than response probability, the density of T regs in lymph nodes was the best biomarker with a highest RIS of 0.35 (Fig. 5E). The virtual patient subgroup identified based on this biomarker incorporated 79% of responders and 43% of nonresponders from the entire cohort ( fig. S7B). The best 30 biomarkers chosen based on RIS included the fraction of TILs that attained a RIS of 0.24. The fraction of TILs is predicted to be a better biomarker in maximizing RIS compared to PD-L1 in tumor and neo-antigen-specific T cell clones, although the latter biomarkers achieved higher response probabilities. Several biomarkers identified as the best 30 based on response probability were also among the best 30 when ranked on the basis of the RIS (biomarkers highlighted in data S2). This includes biomarkers such as the fraction of cytotoxic T cells in lymph nodes/central compartment, fraction of exhausted cytotoxic T cells in the tumor, and the initial tumor diameter. However, the threshold biomarker levels for the virtual patient subgroup that attains the highest response probability or RIS varied (data S2). Together, individual biomarkers tested were able to achieve a maximum RIS of 0.35 and a maximum response probability of 0.86.

Predictive power of combinations of two biomarkers
As the predictive power of individual biomarkers could be limited, we tested whether any combination of two biomarkers achieved higher predictive power. We generated all possible combinations of two biomarker candidates and generated subgroups of virtual patients on the basis of the levels of both biomarker candidates in the combination. The combination of the fraction of cytotoxic T cells in the central compartment with the density of APCs in lymph node achieved a response probability of 1 (Fig. 6A). Similarly, the fraction of cytotoxic T cells in lymph nodes combined with the density of APCs in lymph node also attained a response probability of 1 (Fig. 6A). When the biomarker threshold levels for virtual subgroups were chosen to maximize response probability, the fraction of responders from the entire cohort included in the subgroups were low as observed for single biomarkers ( fig. S8A).
When biomarker combinations were ranked on the basis of the RIS, the highest RIS attained was 0.43 for the combination of immune cell fraction in metastatic tumors and the density of T regs in lymph nodes (Fig. 6B). On the basis of the levels of these biomarker combinations, 72% of responders in the overall cohort were included in the patient subgroup, but only 29% of nonresponders were included ( fig. S8B). For all the biomarker combinations, we calculated the percentage change in RIS/response probability in combination compared to the best individual biomarker in that combination. The increase in predictive power varied greatly among pairs of biomarkers but was considerably enhanced compared to individual ones, suggesting a higher predictive power of biomarker combinations ( fig. S9). Despite only small differences in the response probability/RIS attained among top-ranked biomarkers, the differences were statistically significant (figs. S10 and S11).

Dependence of response rate on model parameters
To identify the most influential parameters in the QSP model, we chose multiple cutoffs for each parameter in the allowable range and generated virtual patient subgroups with parameter values above and below chosen cutoffs. Similar to the biomarker analysis discussed in the previous section, parameters were ranked on the basis of the maximum response probability attained in the subgroups generated. Maximum response probability attained can be considered as a measure of how much the response probability is increased by a particular parameter in a patient subset compared to the entire patient cohort. Therefore, parameters attaining highest response probability were considered as the most influential in determining response rates. Killing rate constant of cancer cells by cytotoxic T cells was the most influential parameter followed by parameters such as macrophage recruitment rate to the other metastatic tumor, tumor vasculature growth rate, and dissociation constant of peptide-major histocompatibility complex (MHC) (Fig. 7A). While increasing the threshold levels for parameters such as cancer cell killing rate, the macrophage recruitment rate resulted in increased fraction of responders; the tumor vasculature growth rate showed the opposite trend (Fig. 7, B to E). We chose 18 model parameters that were fixed among virtual patients and tested the influence on the overall response rates by perturbing the parameter values (50% decrease to 50% increase) ( fig. S12). Among the tested parameters, overall response rates were sensitive to parameters including half-maximal CCL2 level for MDSC recruitment, T cell diversity, rate constant for the degradation of CCL2, rate constant of dead cell clearance from the tumor compartments, rate constant of T cell death by T regs , and the concentration of self-antigen in cancer cells. Response statuses of a notable fraction of patients were shifted from progressive disease (PD) to stable disease (SD) or partial response (PR) by increasing parameters such as T cell diversity, half-maximal CCL2 level for MDSC recruitment, and rate constant for the degradation of CCL2 ( fig. S13). This suggests that a lower T cell diversity or a higher concentration of CCL2 that induces the recruitment of MDSCs and macrophages could be an important factor impeding patient responses. On the other hand, a reverse trend with predominant shift from SD/PR to PD was observed by increasing the values of parameters such as the rate constant of T cell death by T regs and the concentration of self-antigen in cancer cells.
Heterogeneity in the response of metastatic tumors to pembrolizumab therapy Next, we considered a different cohort in which virtual patients are assumed to have intact primary tumor in addition to metastatic tumors. We sought to explore whether the differences in immune microenvironment of primary and metastatic tumors assumed in this study are sufficient to confer differences in response to pembrolizumab therapy. We simulated the response of virtual patients to pembrolizumab monotherapy by administering 200 mg of pembrolizumab every 3 weeks. We compared the increase in primary and metastatic tumor load 400 days after the start of pembrolizumab treatment (Fig. 8). The median primary tumor response to pembrolizumab was higher compared to lung and other metastatic tumors (Fig. 8, A and C). Lung metastatic tumors also showed a significantly better response compared to other metastatic tumors (Fig. 8, A and C). We also observed a large heterogeneity in the response of primary and metastatic tumors (Fig. 8B). Notably, the heterogeneity was higher in both lung and other metastatic tumors compared to primary tumors (Fig. 8B).

DISCUSSION
The treatment of TNBC is challenging owing to its high propensity of metastatic spread and the establishment of metastatic tumors at distant sites such as the lung, liver, brain, and bones (7,8,10). The challenge of developing effective therapeutic strategies is further exacerbated by the intraindividual heterogeneity of metastatic tumors that ultimately leads to heterogeneity in treatment responses. Pembrolizumab monotherapy has demonstrated durable responses in small numbers of patients with previously treated metastatic TNBC (22). However, biomarkers that reliably identify the subgroup of patients who are most likely to achieve therapeutic benefit are now lacking.
In this study, we developed a QSP model for metastatic TNBC that incorporates multiple metastatic tumors to address these challenges in investigating the effects of pembrolizumab. Because lung metastases are observed in most patients with TNBC, we explicitly incorporated metastatic lung tumors into the model. Typically, QSP model development and parameterization rely on data from clinical, preclinical, and in vitro studies. Given the sparsity of such data for metastatic tumors from specific sites, we used transcriptomic data as suggested previously by Lazarou et al. (43) to calibrate the differences in immune cell content of the lung and other metastatic tumor site relative to the primary breast tumor. A virtual clinical trial simulation with the calibrated model was able to capture heterogeneity in the response of tumors within individual patients.
Simulations predicted that pembrolizumab monotherapy not only increased the tumor infiltration of immune-activating cells but also increased immune-suppressive cell types and cytokines. These predictions suggest that a combination of pembrolizumab with therapies targeting immune-suppressive cell types/cytokines such as M2 macrophages, MDSCs, TGF-β, or IL-10 could potentially be beneficial in further increasing overall response rates. We observed that pembrolizumab treatment had similar qualitative effects on most cell types and molecular species in responders and nonresponders except M1 macrophages and IL-12 levels.
Higher levels of biomarkers such as PD-L1 expression, TILs, and tumor mutational burden have been associated with better response (44). We observed an increase in response rates, when patient subgroups with increasing thresholds of PD-L1 expression were analyzed, a trend that was also seen in the KEYNOTE-119 trial (22,27). While PD-L1 expression is a commonly used biomarker in metastatic TNBC, there is a need for better biomarkers because a substantial fraction of patients with high PD-L1 expression are nonresponders and some patients with low PD-L1 expression may respond to pembrolizumab monotherapy (45). We predicted several biomarkers that might potentially confer higher predictive power, such as the richness of cancer clones, the fraction of cytotoxic T cells in lymph node and blood, the fraction of exhausted cytotoxic T cells in metastatic tumors, and the diversity and evenness of cytotoxic T cells in lymph nodes. These predictions suggest that, when virtual patients are selected on the basis of the identified biomarkers with an appropriate threshold biomarker level, the highest response probability that could be achieved in the patient subgroup is 0.86. Because of limited data availability on biomarker analysis, a rigorous quantitative validation of the model could not be performed, and our work generates hypotheses that need to be clinically validated. The proposed predictive biomarkers include those from biopsies and blood samples and can be measured; while biomarkers such as the densities of immune cells can be measured by flow cytometry with appropriate markers, biomarkers such as the richness of cancer clones and the diversity or evenness of cytotoxic T cells

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E
would require genomic and T cell receptor (TCR) sequencing. The presence of predictive biomarkers from blood samples among toprated biomarkers suggests the feasibility of testing biomarkers in cases where obtaining biopsies might be difficult depending on the site of metastases.
We found that biomarkers that resulted in higher response rates in a chosen subset of patients did not necessarily include a large proportion of responders from the population. In addition, only a small fraction of responders from the entire cohort were selected when biomarkers identified to solely maximize response probability were used to select virtual patient subgroups. We thus introduced a metric called RIS to rank biomarkers based on their ability to stratify the responders and nonresponders in the entire patient cohort. While the number of responders was higher in virtual patient subsets selected based on RIS, even a combination of two biomarkers was not able to effectively eliminate the inclusion of nonresponders. Considering the complexity of the tumor microenvironment, we propose that the combination of more than two biomarkers might be required to effectively stratify responders and nonresponders and that future studies focusing on identifying the minimal set of biomarkers required would be invaluable.
Parameters of the QSP model presented were estimated using available in vitro, in vivo, and clinical data. While 12% of the parameters do not have strong experimental support, these parameters were extensively calibrated to fit the clinical response data. However, considering the differences in experimental conditions and tumor samples used, estimated parameters should only be treated as rough estimates that could be refined upon the availability of more data. Sensitivity analysis has identified parameters that, when perturbed, could shift the response status of a substantial fraction of virtual patients. Thus, note that the predictions of the model could rely on parameter values and that model predictions need confirmation from clinical studies. Availability of more data for model calibration and validation would further solidify the ability of the model to generate prospective predictions. Among the model parameters, a subset (127 parameters) was chosen to vary between virtual patients to account for the heterogeneity in patient responses and the uncertainty in parameter values. Compared to previous studies (33,36), a greater number of parameters were varied in this study to account for tumor-tumor heterogeneity within a patient, in addition to interindividual variability. This includes model parameters known to be highly variable among patients such as the initial tumor diameter, PD-L1 expression, and tumor growth rate constant. Parameter distributions for sampling these parameters were estimated by fitting the clinical response data to ensure that there is adequate variability in virtual patient responses. The choice of parameters to vary and their ranges/distributions need to be fine-tuned upon the availability of patient-level data to ensure that the virtual patients generated are realistic.
In this study, we focused on predicting the responses in patients with previously treated metastatic TNBC. The combination of chemotherapy with pembrolizumab was tested in the first-line setting and is approved for the treatment of patients with metastatic PD-L1 + TNBC (25,26). The model and workflow presented in this study could be used to simulate the combination of chemotherapy with pembrolizumab upon the addition of modules for pharmacokinetics and pharmacodynamics of chemotherapies and recalibration of the model parameters to represent the first-line setting. Considering the higher incidence of adverse effects when chemotherapy is combined with pembrolizumab, careful selection of patients for either pembrolizumab monotherapy or pembrolizumab in combination with chemotherapy would also be beneficial in the first-line setting. Similarly, complete disappearance of metastatic tumors was observed in a subset of patients with TNBC treated with atezolizumab alone (46,47) or in combination with nab-paclitaxel (48), suggesting that the identification of appropriate biomarkers would be highly beneficial for patients treated with these drugs (where available) as well. Thus, a comparative in silico analysis of chemo-immunotherapy strategies in the same virtual patient cohort would suggest whether monotherapy with pembrolizumab or atezolizumab could be beneficial for individual patients, or whether they would derive clinical benefit only if chemotherapy was combined with immunotherapy.
A previous study (49) demonstrated that lung metastatic tumors have higher immunogenicity compared to metastatic tumors in the brain, liver, and bones, suggesting that lung metastases might be ideal targets for immunotherapy strategies. However, not all metastatic tumors in the lung show high immunogenicity, suggesting a heterogeneity in treatment responses (49). Consistently, we observed that lung metastatic tumors had a large heterogeneity and, when compared to other metastatic tumors, had a better median response to pembrolizumab. Our simulation results also predicted a relatively resistant nature of metastatic tumors compared to primary tumors that could be a hindering factor in the successful treatment of metastatic TNBC with immunotherapy. It is well known that cold tumors might be relatively resistant to immunotherapy compared to hot tumors, and, in such cases, other treatment options such as chemo-immunotherapy are recommended (50,51).
One of the limitations of our study is that only target lesions were simulated and nontarget/new lesions were neglected in the response evaluation. In clinical trials, the appearance of new metastatic lesions or unequivocal progression of nontarget lesions is considered PD. In the QSP model, we did not consider the detailed process of metastasis and directly introduced cancer cells into tumor compartments at specified time points that were randomly sampled from distributions established using clinical data. This simplified approach minimizes the number of parameters required. However, upon availability of sufficient data to calibrate the model parameters, the detailed process of metastasis could be incorporated to enable a mechanistic understanding of the effect of pembrolizumab on the appearance of new metastatic lesions. Furthermore, explicitly incorporating metastatic tumors other than lung metastases could be of interest in the future. For instance, liver metastases have been associated with systemic suppression of the immune response and poor responses to immunotherapy (52,53). Investigating the response of liver metastatic tumors to various therapeutic strategies would be of importance in exploring treatment options for patients with liver metastases. The availability of patient-level data from clinical trials would facilitate such model advancements. Interaction between tumors due to cytotoxic T cell trafficking was observed to be minimal in the majority of patients in this study (fig. S14). The QSP model developed here can be extended to include mechanisms mediating tumor-tumor interactions for the investigation of abscopal effects seen in the context of local therapy.
Immune cell types such as cytotoxic T cells, helper T cells, T regs , APCs, MDSCs, and macrophages that have been relatively well studied in the TNBC microenvironment were included in the model. Apart from the cell types incorporated in the model, estimates of xCell (54), a gene signature-based cell composition inference tool predicted that lung metastatic tumors tend to have higher stromal cell content compared to the primary tumor ( fig. S1D).
Cancer-associated fibroblasts (CAFs) are a highly heterogenous, prevalent stromal cell type in the tumor microenvironment known to promote angiogenesis, modulate chemoresistance, and facilitate metastasis (55,56). A systemic analysis of stromal cellderived cytokines in TNBC has also identified numerous potential therapeutic targets (57). Apart from the stromal cells, prior analysis of clinical samples has suggested the presence of tertiary lymphoid structures (TLSs) with or without active germinal centers in TNBC tumors (15). The presence of mature TLS with germinal centers has been associated with prognosis in several cancer types (58)(59)(60). These observations suggest that cell types such as CAFs and B cells might affect responses to immunotherapy. Subpopulations of cell types, such as T regs and MDSCs, with distinct phenotypic characteristics and gene expression, have been identified and potentially have different functional capabilities (61,62). Incorporation of these cell types/subtypes in the QSP model in future work would enable testing whether these cell types could potentially act as biomarkers in response prediction to pembrolizumab or other therapies.
In this study, we assumed that primary and metastatic tumors differ mainly in immune cell composition. Despite this assumption, heterogeneity between individual tumors in characteristics other than immune cell composition was considered by randomly sampling parameter values such as clonal diversity and neo-antigen burden for every tumor in a virtual patient. Identical parameter distributions were used for primary and metastatic tumors. Comparing the differences between primary and metastatic tumors in clonal composition and neo-antigens using appropriate high-throughput data from clinical samples of patients with TNBC could help better inform the QSP model developed in this study. While interindividual variability and intraindividual heterogeneity in terms of tumortumor heterogeneity are accounted for in this study, our QSP model ignores intratumoral heterogeneity and assumes that tumors are well-mixed. Spatial immune architecture such as the proximity of certain cell types has been shown to be associated with disease prognosis (63). Agent-based models have been developed and integrated with QSP models to explore the effects of intratumoral heterogeneity on treatment responses (64)(65)(66). Nevertheless, the QSP platform developed in this study could be integrated with such spatial models or adapted for other metastatic cancers to investigate existing and emerging treatment strategies in silico.

QSP model of metastatic TNBC
We developed a QSP model for metastatic TNBC using components of a previously developed QSP model with a single TNBC tumor (33). The extended QSP model consists of metastatic tumor compartments and associated lymph node compartments in addition to the central, peripheral, and primary tumor and primary tumor-draining lymph node compartments considered in the previous study (33). Metastatic tumor compartments include two lung metastatic tumor compartments and a single compartment to represent metastatic tumors other than lung metastases, referred to as other metastatic tumor compartment (Fig. 1). Multiple lung metastases were assumed because it is the most common site of metastasis. As per RECIST v1.1 (67), a maximum of two lesions per organ are chosen as target lesions for quantitative assessment. Therefore, we consider two lung metastatic tumor compartments in the QSP model as only target lesions are simulated in this study. Both lung metastatic tumor compartments are connected to a lung-draining lymph node compartment, and other metastatic tumor compartment is connected to an associated other lymph node compartment via lymphatic system. Each lymph node compartment in the QSP model represents the sum total of lymph nodes draining a particular organ. The number of metastatic tumor compartments is not a limitation of the model; it can be readily extended. The model consists of 641 equations and 737 parameters and is implemented using the SimBiology toolbox of MATLAB (MathWorks, Natick, MA). Here, we discuss extensions introduced in this study and provide only a brief description of the components of the base model. Detailed description of model equations describing dynamics of individual species is provided in (33). The full set of equations, chemical reactions, and model parameters is presented in the Supplementary Materials; the MATLAB code is available at http://dx.doi.org/10.17632/ r46rk4vwdv.1 to ensure reproducibility.
As in the previous model, the QSP model with metastases incorporates multiple modules representing the dynamics of a species or a collection of associated species in the tumor microenvironment. These modules represent the dynamics of cell types such as cancer cells, APCs, macrophages, MDSCs, T regs , cytotoxic T cells, and helper T cells; cytokines such as IL-2, IL-10, IL-12, IFN-γ, TGF-β, and CCL2; angiogenic factors; and checkpoint molecules expressed on immune cells and cancer cells. Mechanisms involved in the dynamics of every model species in metastatic tumor compartments and associated lymph node compartments are assumed to be the same as the primary counterparts.

Brief description of the base model
Dynamics of naïve CD4 + and CD8 + T cells are modeled by considering the trafficking of naïve T cells between the central, peripheral, and lymph node compartments (68) and proliferation in the peripheral and lymph node compartments. Few cancer cells are incorporated into the primary tumor compartment at the start of the simulation, while metastatic tumor compartments are seeded at later time points. Dynamics of cancer cells is modeled by a modified Gompertzian growth pattern, and the dynamics of tumor vasculature is modeled as changes in tumor carrying capacity, which is defined as the maximum number of cancer cells that can be supported at a given time point. Cancer cells are assumed to secrete angiogenic factors that increase the tumor carrying capacity. Natural death or apoptosis of cancer cells is modeled as a first-order reaction and results in the release of tumor-associated neo-antigens and selfantigens into the tumor compartment.
Uptake of tumor-derived neo-antigens and self-antigens by APCs and subsequent maturation results in the migration of mature APCs to the tumor-draining lymph node compartment. Cytokine-mediated effects on APCs are modeled such that IL-12 facilitates the maturation of APCs and that IL-10 inhibits APC maturation. Detailed mechanisms of antigen processing and presentation including cleavage of proteins into peptides in the intracellular vesicles, binding of peptides to MHC molecules, and transport to the cell surface are modeled. Activation of naïve T cells depends on the extent of TCR ligation by peptide-MHC on APCs and is implemented as a Hill function. The number of divisions of activated T cells is approximated as a linear sum of contributions from: (i) TCR ligation, (ii) CD28 engagement that is competitively inhibited by CTLA-4, and (iii) IL-2 derived from helper T cells (69). A separate compartment representing synapse between T cell and APC is considered to model the interactions of CD28 and CTLA-4 with their shared ligands CD80 and CD86. Following activation and transport into the central compartment, infiltration of T regs , cytotoxic T cells, and helper T cells into the tumor compartment is assumed to depend on the tumor volume and volume fraction of vascular space in tumor. Tumor-infiltrating cytotoxic T cells kill cancer cells with a rate depending on the ratio of cytotoxic T cells and cancer cells and result in the enhanced release of tumor-associated antigens. Tumor-infiltrating cytotoxic T cells and helper T cells are assumed to become exhausted by the interaction of PD-1 with ligands on cancer cells and by the action of IL-10 from T regs .
Cancer cells secrete CCL2 that is assumed to mediate the recruitment of M1 macrophages and MDSCs into the tumor compartment. M1 macrophages undergo reversible polarization to M2 macrophages. M1 macrophages are the source of IL-12, and M2 macrophages secrete angiogenic factors, TGF-β, and IL-10. Cytokines IL-10 and TGF-β promote the polarization of M1 to M2 macrophages, whereas cytokines IL-12 and IFN-γ promote the polarization of M2 to M1 macrophages. Phagocytosis of cancer cells by M1 macrophages is inhibited by IL-10 and the ligation of immune checkpoints PD-1 and signal regulatory protein α (SIRPα). MDSCs are assumed to release arginase I (Arg-I) and nitric oxide (NO). Cytotoxic activity of T cells is inhibited by TGF-β, Arg-I, and NO. TGF-β and Arg-I facilitate the trans-differentiation of helper T cells to T regs in the tumor.
The QSP model includes pharmacokinetics and pharmacodynamics of the anti-PD-1 antibody pembrolizumab. The antibody is directly administered into the central compartment mimicking intravenous infusion, and the pharmacokinetics incorporates antibody clearance from the central compartment, transport between central and peripheral/tumor compartments, and the transport from tumor to tumor-draining lymph node compartments. Pharmacodynamics of anti-PD-1 antibody is modeled by the binding of the antibody to PD-1 receptor on cytotoxic T cells and macrophages, which blocks the interactions of PD-1 with PD-L1 and PD-L2. Subsequent reduction in the amount of ligand-bound PD-1 decreases the inhibitory action of PD-1 on T cell and macrophagemediated killing/phagocytosis of cancer cells that are both modeled as Hill functions.

Extensions to clonal diversity and neo-epitopes
To account for the clonal heterogeneity of individual tumors, five cancer clones are considered. Apart from the number of cells of each cancer clone, individual cancer clones are assumed to differ in growth rate and neo-epitope expression. The presence of eight neo-epitopes (70) and cytotoxic T cells with different specificities corresponding to each neo-epitope is considered. Assuming that each neo-epitope is recognized by multiple TCR clonotypes, cytotoxic T cell population specific to each neo-epitope is assumed to consist of multiple T cell clonotypes. We also assume that cytotoxic T cells of any neo-epitope specificity can recognize and kill every cancer clone; however, the number of activated cytotoxic T cells specific to each epitope is dependent on the concentration of corresponding peptides presented by mature APCs.

Model calibration
Model parameters were estimated in the previous versions of the QSP model using in vitro, preclinical, and clinical data or directly obtained from the literature. Tumor growth parameters were estimated by fitting the model to tumor growth curve from TNBC xenograft models and were scaled up to humans by allometric scaling. The initial amount of naïve T cells is based on the measurement from healthy individuals, and the rate of thymic export of naïve T cells is estimated using average ages of patients with breast cancer. Recruitment rate constants of immune cells are estimated to fit the cell density measurements from breast cancer samples. Cell surface expression of checkpoint molecules is obtained from flow cytometry measurements. Secretion of soluble factors such as cytokines is calibrated using the cytokine concentration in breast cancer samples. Expression levels of Arg-I and NO are estimated using in vitro data from breast cancer cells. Pharmacokinetic parameters of therapeutic antibody are estimated to fit the clinical measurements of plasma concentration. These parameters were inherited into the extended QSP model and left unchanged.
However, immune cell migration rate constants and cellular densities for the metastatic tumors are calibrated such that the immune cell content of metastatic tumors is consistent with cell abundance estimations from publicly available transcriptomic data as discussed in the "Calibration of the immune cell content of metastases" section. Parameters such as the seeding time of metastatic tumors, the initial number of cancer cells in tumor compartments, the number of antigen-specific T cell clones, the rate of T cell exhaustion by cancer cells, and the initial diameter of metastases were estimated to fit the clinical response data from KEYNOTE-119 trial as explained in the "Virtual clinical trial" section. The list of model parameters and references are provided in data S1.
The Huang dataset included transcriptomes of metastatic tumors from different sites of a TNBC patient and was used for the comparison of immune cell abundance between metastatic tumors in the lung, liver, brain, and bone (see figs. S1A and S2, A and B). After excluding low-purity samples [supplementary text of (11)], this dataset included seven lung, four liver, three brain, and three bone metastatic tumor samples. The Siegel dataset included transcriptomes of 39 metastatic tumors from sites such as the lung, liver, brain, bone, kidney, skin, adrenal and soft tissue, and matched primary tumors from nine patients with TNBC. We used a gene signature-based cell composition inference tool called xCell (54) to calculate immune and stroma score that provide the overall abundance of immune cells and stromal cells, respectively, in each sample using transcriptomic data.

Calibration of the immune cell content of metastases
For the QSP model calibration, only the abundance estimates of EPIC and quanTIseq were used as these estimates can be interpreted as cell fractions/absolute cell counts (76) and can be directly compared with cell fractions from the simulations. The relative abundance of different cell types in metastatic tumors was calculated with respect to the matched primary tumor from the simulation results and compared with the relative abundance estimates from transcriptomic data. All metastatic tumor samples excluding the lung, brain, and lymph nodes were lumped together to estimate the immune cell composition of other metastatic tumors. Model parameters such as the recruitment rate constants of macrophages and MDSCs, the polarization rate constant of M1 to M2 macrophages, and the baseline density of APCs were estimated such that the abundance of cell types from simulations was consistent with estimates from the RNA-seq data. We use the relative abundance of immune cells in metastases with respect to primary tumor as the relative abundance estimates are expected to be less sensitive to the gene expression platform and normalization techniques used (76).
The transcriptomic data used to estimate immune cell abundances were collected from metastatic tumor samples from patients heavily pretreated with different chemotherapies. By calibrating the immune cell composition of metastases before the start of pembrolizumab treatment using these estimates, we assume that metastatic tumors in the model are heavily pretreated with chemotherapies as in the patients enrolled in KEYNOTE-119 trial (22). Thus, the effect of pretreatment with chemotherapies on immune cell composition is considered without the explicit simulation with chemotherapies in the model.

Virtual clinical trial
Virtual patients were generated by Latin hypercube sampling of selected parameters, such as the seeding time of metastatic tumors, the initial tumor diameter, and the growth rate constant of cancer clones, from a chosen distribution (data S1); each parameter set corresponds to a virtual patient. Parameters randomly sampled account for intraindividual variability (tumor-to-tumor heterogeneity) in addition to the interindividual variability considered in previous studies (33,36). For each virtual patient, a few cancer cells were introduced into metastatic tumor compartments at randomly generated seeding time points, and the simulation was continued until a metastatic tumor reached the initial target diameter. Values of model species at the end of the simulation were used as initial conditions for treatment simulations. Among the 1000 virtual patients generated, only the virtual patients (816 of 1000) who met the target tumor diameter were considered for treatment simulations. Parameter sets corresponding to virtual patients who did not attain the target tumor diameter were discarded to eliminate potentially unrealistic virtual patients. After the start of virtual treatment, no metastatic tumors were seeded. This results in different number of tumors in simulated virtual patients despite the presence of three metastatic tumor compartments in all virtual patients. While virtual patients were filtered only on the basis of the ability to develop a tumor in this study, the availability of more clinical data would enable filtering virtual patients based on plausible ranges of species concentrations (78).
The virtual clinical trial simulation was performed by administering pembrolizumab at a dose of 200 mg every 3 weeks following the pembrolizumab arm of the KEYNOTE-119 trial (22). Assuming that most patients underwent primary tumor resection before enrollment in the KEYNOTE-119 trial, the initial number of cancer cells in the primary tumor compartment was set to 0 in all virtual patients.
Response of each virtual patient was categorized as complete response (CR)/PR, SD, or PD based on the RECIST v1.1 criteria (67). As in the KEYNOTE-119 clinical trial, a minimum duration of 24 weeks was considered for SD, and time intervals between subsequent observations for response evaluation were set to 9 weeks. Proportion of patients with CR/PR was considered as overall response rate. Only virtual patients who achieved a CR or PR were included in the calculation of the time to response and duration of response. Duration of response in patients who did not relapse until the end of simulation was set as the duration from response to the end of simulation. The total number of patients used in response rate calculations reported from the KEYNOTE-119 trial also includes nonevaluable patients (22). Thus, we renormalized the results from KEYNOTE-119 by excluding nonevaluable patients to enable direct comparison with the virtual clinical trial. Bootstrap sampling was performed to calculate 95% confidence intervals for response rates, duration of response, and time to response of the virtual patient population. We assumed that the percentage of patients with lung metastases in pembrolizumab arm is 65%, which is the percentage of patients with lung metastases in KEYNOTE-119 trial at baseline (27). Median seeding time of metastatic tumors was estimated to fit the percentage of patients with lung metastases in the virtual clinical trial. Because of the lack of data on the baseline metastatic tumor diameter, this parameter was estimated to fit the response rates resulting in a median tumor diameter of 1.65 cm that is within clinically reported ranges of the sizes of metastatic tumors (11).

Biomarker analysis
We classified the virtual patients as responders and nonresponders based on the response status. In the analysis of treatment effects and biomarkers, patients demonstrating CR/PR or SD were considered as responders, whereas nonresponders included patients with PD. To identify and evaluate biomarkers, we considered two different metrics in this study: (i) response probability defined as the number of responders in a chosen patient subgroup divided by the total number of patients in the subgroup; (ii) RIS is defined as follows

RIS ¼
No. of responders in a subgroup Total no. of responders in the entire cohort À No. of nonresponders in a subgroup Total no. of nonresponders in the entire cohort ð1Þ We selected densities/concentrations of different cellular and molecular species from central, lymph node, and tumor compartments as biomarker candidates. In addition, we also considered derived quantities such as the ratio of M1 and M2 macrophages. Taking advantage of the consideration of multiple cancer clones, neo-antigens, and cytotoxic T cell specificities in the QSP model, we considered diversity indices (79,80) of cancer cells and cytotoxic T cells as biomarker candidates.
Richness (S) is defined as the number of unique cancer clones or unique neo-antigen specificities among cytotoxic T cells. Shannon index (H ) and evenness (J ) are calculated using the following equations: k represents different clones for cancer cells or different neoantigen specificities for cytotoxic T cells; p k represents proportion of cancer cells belonging to clone k or proportion of cytotoxic T cells with neo-antigen specificity k For biomarker candidates from tumor or tumor-draining lymph node compartments, average values from all simulated metastatic tumors or tumor-draining lymph nodes, respectively, were calculated. For each biomarker candidate, the range of baseline level was determined in virtual patients previously generated for treatment simulations, and eight uniformly spaced threshold values/cutoffs were chosen, spanning the entire range of biomarker levels. Virtual patients with biomarker candidate values above and below chosen threshold values were selected as different subgroups. This results in virtual patient subgroups that are not mutually exclusive. Subgroups with total number of patients less than 20 were discarded. For each subgroup of virtual patients, response probability and RIS were calculated. Subgroup with highest response probability and RIS was considered as the best subgroup for each biomarker candidate. Biomarker candidates were ranked separately on the basis of the values of the two metrics, response probability and RIS, to identify potential biomarkers. To identify combinations of two biomarkers, all possible combinations of two biomarker candidates were generated. For each biomarker candidate combination, virtual patient subgroups were chosen on the basis of threshold values of both biomarker candidates. Biomarker combinations were then ranked on the basis of the highest response probability or RIS achieved in the best subgroup identified. Biomarker threshold values in virtual patient subgroups with highest response probability and RIS are provided in data S2.

Supplementary Materials
This PDF file includes: Figs. S1 to S14 Legends for data S1 and S2 Other Supplementary Material for this manuscript includes the following: Data S1 and S2 View/request a protocol for this paper from Bio-protocol.